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We study the ±J random-plaquette Z2 gauge model (RPGM) in three spatial dimensions, a 
three-dimensional analog of the two-dimensional ±J random-bond Ising model (RBIM). The model 
is a pure Z2 gauge theory in which randomly chosen plaquettes (occuring with concentration p) have 
couplings with the "wrong sign" so that magnetic flux is energetically favored on these plaquettes. 
Excitations of the model are one-dimensional "flux tubes" that terminate at "magnetic monopoles" 
located inside lattice cubes that contain an odd number of wrong-sign plaquettes. Electric con- 
finement can be driven by thermal fluctuations of the flux tubes, by the quenched background of 



magnetic monopoles, or by a combination of the two. Like the RBIM, the RPGM has enhanced 
symmetry along a "Nishimori line" in the p-T plane (where T is the temperature) . The critical con- 
centration p c of wrong-sign plaquettes at the confinement-Higgs phase transition along the Nishimori 

(N ■ line can be identified with the accuracy threshold for robust storage of quantum information using 

topological error-correcting codes: if qubit phase errors, qubit bit-flip errors, and errors in the mea- 
surement of local check operators all occur at rates below p c , then encoded quantum information 
can be protected perfectly from damage in the limit of a large code block. Through Monte-Carlo 

\q • simulations, we measure p c o, the critical concentration along the T — axis (a lower bound on p c ), 

I finding p c o — .0293 ± .0002. We also measure the critical concentration of antiferromagnetic bonds 

in the two-dimensional RBIM on the T — axis, finding p c o — .1031 ± .0001. Our value of p c o is 
incompatible with the value of p c = .1093 ± .0002 found in earlier numerical studies of the RBIM, in 
^ , disagreement with the conjecture that the phase boundary of the RBIM is vertical (parallel to the 

T axis) below the Nishimori line. The model can be generalized to a rank-r antisymmetric tensor 
field in d dimensions, in the presence of quenched disorder. 



interactions among qubits arranged in a two-dimensional 
block, and the protected information is associated with 
global topological properties of the quantum state of the 
block. If the error rate is small, then the topological 
properties of the code block are well protected, and er- 
ror recovery succeeds with a probability that rapidly ap- 
proaches one in the limit of a large code block. But if the 
error rate is above a critical value, the accuracy threshold, 
then quantum error correction is ineffective. 

In [9], a precise connection was established between 
the accuracy threshold achievable with surface codes and 
the confinement-Higgs transition in a three-dimensional 
Zi lattice gauge model with quenched randomness. The 
model has two parameters: the temperature T and the 
concentration p of "wrong-sign" plaquettes. On wrong- 
sign plaquettes (which are analogous to antiferromag- 
netic bonds in a spin system) it is energetically favorable 
for the Z2 magnetic flux to be nontrivial. In the mapping 
between quantum error recovery and the gauge model, 
the quenched fluctuations correspond to the actual er- 
rors introduced by the environment; these impose sites 
of frustration, magnetic monopoles, corresponding to an 
"error syndrome" that can be measured by executing a 
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I. INTRODUCTION 



Spin systems with quenched randomness have been ex- 
tensively studied, leading to valuable insights that ap- 
ply to (for example) spin glass materials, quantum Hall 
systems, associative memory, error-correcting codes, and 
combinatorial optimization problems [1-3]. Gauge sys- 
tems with quenched randomness, which have received 
q-i comparatively little attention, will be studied in this pa- 
per. 

The gauge models we consider are intrinsically interest- 
ing because they provide another class of simple systems 
with disorder-driven phase transitions. But our inves- 
tigation of these models has a more specific motivation 
connected to the theory of quantum error correction. 

In practice, coherent quantum states rapidly decohere 
due to uncontrollable interactions with the environment. 
But in principle, if the quantum information is cleverly 
encoded [6,7], it can be stabilized and preserved using 
fault-tolerant recovery protocols [8]. Kitaev [4,5] pro- 
posed a particularly promising class of quantum error- 
correcting codes (surface codes) in which the quantum 
processing required for error recovery involves only local 
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suitable quantum circuit. Thermally fluctuating mag- 
netic flux tubes, which terminate at magnetic monopoles, 
correspond to the ensemble of possible error patterns that 
could generate a particular error syndrome. (The tem- 
perature T is tied to the strength p of the quenched fluc- 
tuations through a Nishimori relation [10].) When the 
disorder is weak and the temperature low (correspond- 
ing to a small error rate), the system is in a magnetically 
ordered Higgs phase. In the surface code, magnetic or- 
der means that all likely error patterns that might have 
produced the observed error syndrome are topologically 
equivalent, so that the topologically encoded information 
resists damage. But at a critical value p c of the disorder 
strength (and a temperature determined by Nishimori's 
relation), magnetic flux tubes condense and the system 
enters the magnetically disordered confinement phase. In 
the surface code, magnetic disorder means that the error 
syndrome cannot point to likely error patterns belong- 
ing to a unique topological class; therefore topologically 
encoded information is vulnerable to damage. 

Although the code block is two dimensional, the gauge 
model is three dimensional because one dimension repre- 
sents time. Time enters the analysis of recovery because 
measurements of the error syndrome might themselves 
be faulty; therefore measurements must be repeated on 
many successive time slices if they are to provide reli- 
able information about the errors that have afflicted the 
code block. If qubit phase errors, qubit bit-flip errors, 
and errors in the measurement of local check operators 
all occur at rates below p c , then encoded quantum infor- 
mation can be protected perfectly from damage in the 
limit of a large code block. As we consider more and 
more reliable measurements of the syndrome, the corre- 
sponding three-dimensional gauge model becomes more 
and more anisotropic, reducing in the limit of perfect 
measurements to the two-dimensional random-bond Ising 
model. 

The numerical value p c of the accuracy threshold is of 
considerable interest, since it characterizes how reliably 
quantum hardware must perform in order for a quantum 
memory to be robust. In the three-dimensional Z2 gauge 
model, p c is the value of the wrong-sign plaquette concen- 
tration where the confmcment-Higgs boundary crosses 
the Nishimori line in the p-T plane. A lower bound on 
p c is provided by the critical concentration p c0 on the 
T = axis. In [9], an analytic argument established that 
Pco > .0114. In this paper we report on a numerical 
calculation that finds p c o = .0293 ± .0002. 

In the case where the error syndrome can be measured 
flawlessly, the critical error rate is given by the critical 
antifcrromagnetic bond concentration on the Nishimori 
line of the two-dimensional random-bond Ising model 
(RBIM). Numerical calculations performed earlier by 
other authors [11,12] have established p c = .1093±.0002. 
According to a conjecture of Nishimori [13] and Kitatani 
[14], this value of p c should agree with the critical bond 
concentration p c o of the 2D RBIM on the T = axis. 
The same reasoning that motivates this conjecture for 



the RBIM indicates that p c = p c o for the 3D random- 
plaquette gauge model (RPGM) as well. However, we 
have calculated p c o in the 2D RBIM numerically, finding 
Pco = .1031 ± .0001. Our value of p c0 agrees with an 
earlier numerical calculation by Kawashima and Ricgcr 
[23], but disagrees with the conjecture that p c = p c o- 

In Sec. II we describe in more detail the properties 
of the 2D RBIM and the 3D RPGM, emphasizing the 
importance of the Nishimori line and the inferences that 
can be made about the behavior of order parameters on 
this line. Section III reviews the connection between the 
models and error recovery using surface codes. Our nu- 
merical results for p C Q and for the critical exponent vq at 
the T = critical point are presented in Sec. IV. Section 
V summarizes our conclusions. 



II. MODELS 

A. Random-bond Ising model 

The two-dimensional ±J random-bond Ising model 
(RBIM) has a much studied multicritical point at which 
both the temperature and the strength of quenched dis- 
order are nonzero. This model is an Ising spin system on 
a square lattice, with a variable Sj = ±1 residing at each 
lattice site i. Its Hamiltonian is 

H = -J^njSiSj , (l) 

where J is the strength of the coupling between neighbor- 
ing spins, and Tjj = ±1 is a quenched random variable. 
(That is, nj depends on what sample of the system is 
selected from a certain ensemble, but is not subject to 
thermal fluctuations.) The r^-'s are independently and 
identically distributed, with the antiferromagnetic choice 
Tij = — 1 (favoring that neighboring spins antialign) oc- 
curing with probability p, and the ferromagnetic choice 
Tij=+1 (favoring that neighboring spins align) occuring 
with probability 1 —p. We refer to p as the concentration 
of antiferromagnetic bonds, or simply the bond concen- 
tration. 

The free energy F of the model at inverse temperature 
(3, averaged over samples, is 

[pF(K, r)] Kp =-J2 P(K P , r) In Z(K, r) (2) 

r 

where 




is the partition function for sample r (with K = /3J), 
and 
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P(K p , t) = (2 cosh K p )~ Nb x cxp K p ^ Tij (4) 

V («> / 

is the probability of the sample r; here 



and iV B is the number of bonds. 

The partition function Z(K,t) is invariant under the 
change of variable 

Si ► (TiSi , Tij > <Ti<JjTij , (6) 

where Uj = ±1. Thus r itself has no invariant meaning 
— samples t and r' that differ by the change of variable 
have equivalent physics. The only invariant property of 
t that cannot be modified by such a change of variable 
is the distribution of frustration that r determines. If an 
odd number of the bonds contained in a specified plaque- 
tte have r = — 1 then that plaquette is frustrated — an 
Ising vortex resides at the plaquette. For purposes of vi- 
sualization, we sometimes will find it convenient to define 
the spin model on the dual lattice so that the spins reside 
on plaquettes and the Ising vortices reside on sites. Then 
excited bonds with TijSiSj = — 1 form one-dimensional 
chains that terminate at the frustrated sites. 

Changes of variable define an equivalence relation on 
the set of 2 Nb t configurations: there are the 2 Ns ele- 
ments of each equivalence class (the number of changes of 
variable, where N$ is the number of sites) and there are 
2^ classes (the number of configurations for the Ising 
vortices — note that Nb = 2N$ for a square lattice 
on the 2-torus, and that the number of plaquettes is 
Np = Ns). Denote a distribution of Ising vortices, or 
equivalently an equivalence class of r's, by r\. The prob- 
ability P(K pi 77) of 77 is found by summing P(K p , r) over 
all the representatives of the class; hence 

(2 cosh K p ) NB P(K p ,r]) = (2coshK p ) NB Y,P(K P ,T) 

r£r) 

= CXP \ K P T ij a i a 3 = Z ( K P> V) ■ ( 7 ) 

Apart from a normalization factor, the probability of a 
specified distribution of frustration is given by the par- 
tition function of the model, but with K = [3 J replaced 
by K p . 

In this model, we can define an order parameter 
that distinguishes the ferromagnetic and paramagnetic 
phases. Let 

m 2 (K,K p )= lim [(S i S j ) K } K , (8) 

|»— j|— »oo p 

where (-)k denotes the average over thermal fluctuations, 
[•]k p denotes the average over samples, and denotes 



the distance between site i and site j; then in the ferro- 
magnetic phase m 2 > and in the paramagnetic phase 
to 2 = 0. But the two-point correlation function (SiSj)K 
is not invariant under the change of variable eq. (6), so 
how should to 2 be interpreted? 

Following [9], denote by E the set of bonds that are 
antiferromagnetic (t^ = —1), denote by E' the set of 
excited bonds with r^SiSj = — 1, and denote by D the 
set of bonds with SiSj = — 1 (those such that the neigh- 
boring spins antialign) — see Fig. 1. Then D = E + E' 
is the disjoint union of E and E 1 (containing bonds in E 
or E' but not both). Furthermore, D contains an even 
number of the bonds that meet at any given site; that 
is, D is a cycle, a chain of bonds that has no boundary 
points. The quantity SiSj just measures whether a line 
connecting i and j crosses D an even number {SiSj = 1) 
or an odd number (SiSj = —1) of times. 




-1 -1 -1 I +1 




FIG. 1. The chain E of antiferromagnetic bonds (darkly 
shaded) and the chain E' of excited bonds (lightly shaded) , 
in the two-dimensional random-bond Ising model. Ising spins 
taking values in {±1} reside on plaquettes; Ising vortices 
(boundary points of E) are located on the sites marked by 
filled circles. The bonds of E' comprise a one-dimensional de- 
fect that connects the vortices. The cycle D — E + E' encloses 
a domain of spins with the value —1. 

Now D consists of disjoint "domain walls" that form 
closed loops. If loops that are arbitrarily large appear 
with appreciable weight in the thermal ensemble, then 
the two-point function {SiSj)K decays like exp(— \i— j\/f;) 
— fluctuations far from the sites i and j contribute to the 
correlation function. Thus the spins are disordered and 
m 2 = 0. But if large loops occur only with negligible 
probability, then only fluctuations localized near i and j 
contribute significantly; the spin correlation persists at 
large distances and m 2 > 0. Thus, the order parameter 
probes whether the chain E' of excited bonds can wan- 
der far from the chain E of ferromagnetic bonds; that is, 
whether D = E + E' contains arbitrarily large connected 
closed loops, for typical thermal fluctuations and typical 
samples. 

Nishimori [10] observed that the model has enhanced 
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symmetry properties along a line in the p-T plane (the 
Nishimori line) denned by K = K p or exp(— 2/3 J) = 
p/(l — p). In this case, the antiferromagnetic bond 
chain E and the excited bond chain E' are generated 
by sampling the same probability distribution, subject 
to the constraint that both chains have the same bound- 
ary points. This feature is preserved by renormal- 
ization group transformations, so that renormalization 
group flow preserves Nishimori's line [15]. The Nishi- 
mori point (p c ,T c ) where the Nishimori line crosses the 
ferromagnetic-paramagnetic phase boundary, is a renor- 
malization group fixed point, the model's multicritical 
point. 

When the temperature T is above the Nishimori line, 
excited bonds have a higher concentration than antiferro- 
magnetic bonds, so we may say that thermal fluctuations 
play a more important role than quenched randomness 
in disordering the spins. When T is below the Nishimori 
line, antiferromagnetic bonds are more common than ex- 
cited bonds, and the quenched randomness dominates 
over thermal fluctuations. Right on the Nishimori line, 
the effects of thermal fluctuations and quenched random- 
ness are in balance [16]. 

By invoking the change of variable eq. (6), various 
properties of the model on the Nishimori line can be de- 
rived [3,10]. For example, the internal energy density (or 
"average bond" ) can be computed analytically, 



[n j (S i S j ) Kp } Kp = l-2p 



(9) 



where i and j are neighboring sites; averaged over ther- 
mal fluctuations and samples, the concentration of ex- 
cited bonds is p as one would expect (and the internal 
energy has no singularity at the Nishimori point). Fur- 
thermore, after averaging over disorder, the (2m — l)st 
power of the k-spm correlator has the same value as the 
(2m)th power, for any positive integer m: 



\2m-l 



K„ 



{Si 1 Si. 



c \2m 



K„ 



(10) 



It follows in particular that the spin-glass order parame- 
ter 



q 2 (K p ,K p ) 



lim 



{SiSj) 2 K 



(11) 



coincides with the ferromagnetic order parameter 
m 2 (K p , K p ) along the Nishimori line, reflecting the prop- 
erty that thermal fluctuations and quenched randomness 
have equal strength on this line. 
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FIG. 2. The phase diagram of the random-bond Ising 
model (shown schematically), with the temperature T on 
the vertical axis and the concentration p of antiferromagnetic 
bonds on the horizontal axis. The solid line is the boundary 
between the ferromagnetic (ordered) phase and the param- 
agnetic (disordered) phase. The dotted line is the Nishimori 
line e~ 2,3J = p/(l — p), which crosses the phase boundary at 
the Nishimori point (the heavy black dot). It has been con- 
jectured, but not proven, that the phase boundary from the 
Nishimori point to the p-axis is vertical, as in (a). The nu- 
merics reported in Sec. IV favor the reentrant phase diagram 
shown in (b). The deviation of the critical bond concentration 
p c on the Nishimori line from the critical bond concentration 
PcO on the T — axis has been exaggerated in (b) for clarity. 

Comparing eq. (2) and (7), we see that for K = K p 
the free energy of the model coincides with the Shan- 
non entropy of the distribution of vortices, apart from a 
nonsingular additive term: 

[/3F(K p ,T)} Kp 

= - J2p(K p ,v)^P(K p ,t 1 ) - N B \n (2 cosh K p ) . (12) 

Since the free energy is singular at the Nishimori point 
{p c ,T c ), it follows that the Shannon entropy of frustra- 
tion (which does not depend on the temperature) is sin- 
gular at p — p c [13]. This property led Nishimori to sug- 
gest that the boundary between the ferromagnetic and 
paramagnetic phases occurs at p — p c at sufficiently low 
temperature, and thus that the phase boundary is ver- 
tical in the p-T plane below the Nishimori point, as in 
Fig. 2a. Later, Kitatani [14] arrived at the same conclu- 
sion by a different route, showing that the verticality of 
the phase boundary follows from an "appropriate condi- 
tion." These arguments, while suggestive, do not seem 
compelling to us. There is no known rigorous justifica- 
tion for Kitatani's condition, and no rigorous reason why 
the ferro-para boundary must coincide with the singular- 
ity in the entropy of frustration, even at low tempera- 
ture. Hence we regard the issue of the verticality of the 
phase boundary as still unsettled. Nishimori did argue 
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convincingly that the phase boundary cannot extend to 
any value of p greater than p c [10], and Le Doussal and 
Harris argued that the tangent to the phase boundary 
is vertical at the Nishimori point [15], but these results 
leave open the possibility of a "reentrant" boundary that 
slopes back toward the T axis below the Nishimori point, 
as in Fig. 2b. 

The RBIM can also be defined in d dimensions. Much 
of the above discussion still applies, with minor modifi- 
cations. Consider, for example, d — 3. On the dual lat- 
tice, spins reside on lattice cubes and the bonds become 
plaquettes shared by two neighboring cubes. The set of 
antiferromagnetic bonds E is dual to a two-dimensional 
surface, and its boundary dE consists of one-dimensional 
loops — the Ising strings where the spins are frustrated. 
The set of excited bonds E' is dual to another two- 
dimensional surface that is also bounded by the Ising 
strings: dE' = dE. The spins are disordered if the 
two-cycle D = E + E' contains arbitrarily large closed 
connected surfaces for typical thermal fluctuations and 
typical samples. Similarly, in d dimensions, frustration is 
localized on closed surfaces of dimension d — 2, and the 
thermally fluctuating defects are dimension- (d — 1) sur- 
faces that terminate on the locus of frustration. For any 
d, the model has enhanced symmetry along the Nishimori 
line K = K p , where antiferromagnetic bonds and excited 
bonds are drawn from the same probability distribution. 

In the absence of quenched disorder, the two- 
dimensional Ising model is mapped to itself by a duality 
relation that can be used to infer properties of the critical 
theory. When quenched disorder is introduced, however, 
the two-dimensional random bond Ising model is mapped 
under duality to a model with Boltzmann weights that 
are not positive definite [17], so that it is not easy to draw 
any firm conclusions. 



B. Random-plaquette gauge model 

In the d-dimensional RBIM, excitations have codimcn- 
sion 1 and terminate on a closed surface of codimcn- 
sion 2. The Z 2 random-plaquette gauge model (RPGM) 
is defined in an entirely analogous manner, except that 
the excitations are codimension-2 objects ("magnetic flux 
tubes") that terminate on codimcnsion-3 objects ("mag- 
netic monopoles"). 

More concretely, the variables of the model are Ut = 
±1 residing on each link I of the lattice, and the Hamil- 
tonian is 



H = -jJ2rpU P , 



where J is the coupling strength, 

Up = J[ U t 

eeP 



(13) 



(14) 



is the ^-valued "magnetic flux" through the plaqucttc 
P, and Tp = ±1 is a quenched random variable. The Tp's 
are independently and identically distributed, with the 
"wrong-sign" choice Tp = — 1 (favoring nontrivial flux) 
occuring with probability p, and the "right-sign" choice 
rp=+l (favoring trivial flux) occuring with probability 
1 — p. We refer to p as the concentration of wrong-sign 
plaquettes, or simply the plaquette concentration. 

The free energy F of the model at inverse temperature 
0, averaged over samples, is 

[0F(K, T)] Kp =-J2 P{K P , t) In Z(K, r) (15) 



where 



Z(K, t)=J2 exp (k J2 t pUp^ 



(16) 



is the partition function for sample r (with K = J), 
and 

P(K P , t) = {2coshK p )- Np x exp ^K p ^ ( 17 ) 



is the probability of the sample r; here 



1-p 



(18) 



and Np is the number of plaquettes. 

The partition function Z(K,t) is invariant under the 
change of variable 



t p — > apTp , 



(19) 



where erg = ±1 and op = Y\ eeP ere- While r itself has no 
invariant meaning, r determines a distribution of frus- 
tration that cannot be altered by a change of variable. 
If an odd number of the plaquettes contained in a speci- 
fied cube have r = — 1 then that cube is frustrated — a 
Z2 magnetic monopole resides in the cube. For purposes 
of visualization, we will sometimes find it convenient to 
define the gauge model on the dual lattice so that the 
gauge variables Ut reside on plaquettes, the magnetic flux 
on bonds, and the magnetic monopoles on sites. Then 
excited bonds with TpUp = —1 form one-dimensional 
strings that terminate at monopoles. 

We can define an order parameter that distinguishes 
the Higgs (magnetically ordered) phase and the confine- 
ment (magnetically disordered) phase. Consider the Wil- 
son loop operator associated with a closed loop C (on the 
original lattice, not the dual lattice): 



W{C) 



tec 



(20) 



and consider the behavior of the expectation value of 
W(C), averaged over thermal fluctuations and over sam- 
ples. In the Higgs phase, for a large loop C the Wilson 
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loop operator decays exponentially with the perimeter of 
the loop, 



{(W(C)) 



K 



exp [— [i ■ Perimeter(C)] , (21) 



while in the confinement phase it decays exponentially 
with the area of the minimal surface bounded by C, 



l(W(C)) K ] K ~expl-K-ATe a (C)} 



(22) 



The interpretation is that on the dual lattice the wrong- 
sign plaquettes correspond to a one-chain E bounded by 
magnetic monopoles, and the excited plaquettes corre- 
spond to another one-chain E' with the same boundary; 
hence D = E + E' is a cycle, a sum of disjoint closed "flux 
tubes." If arbitrarily large loops of flux appear with ap- 
preciable weight in the thermal ensemble for typical sam- 
ples, then magnetic fluctuations spanning the entire sur- 
face bounded by C contribute to the expectation value 
of W(C), and the area-law decay results. If large flux 
tubes are suppressed, then only the fluctuations local- 
ized near the loop are important, and the perimeter-law 
decay applies. Thus, the Wilson-loop order parameter 
probes whether the chain E' of excited plaquettes can 
wander far from the chain E of wrong-sign plaquettes; 
that is, whether D = E + E' contains arbitrarily large 
connected closed loops. 

The one-chain E bounded by the magnetic monopoles 
is analogous to a ^-valued Dirac string — the change of 
variable eq. (19) deforms the strings while leaving invari- 
ant the boundary of E (the locations of the monopoles) . 
One should notice that these strings are not invisible to 
our Wilson loop operator; that is W(C) is not invariant 
under the change of variable. It is possible to modify 
W(C) to obtain an invariant object [18], but that would 
not be appropriate if the order parameter is supposed to 
probe the extent to which the thermally fluctuating de- 
fects (the excited plaquettes) depart from the quenched 
disorder (the Dirac strings). 

Like the RBIM, the RPGM has enhanced symmetry 
on the Nishimori line K — K p , and the change of vari- 
able eq. (19) may be invoked to derive properties of the 
model on this line. The Nishimori line is preserved by 
rcnormalization group flow, and crosses the confinement - 
Higgs boundary at a multicritical point (p c ,T c ). The 
internal energy (or average plaquette) can be computed 
on this line, 



[r P (U P ) Kp ] =l-2p 



(23) 



(excited plaquettes have concentration p) and for each 
positive integer m, the (2m — l)'st power of W(C) and 
the 2m'th power are equal when averaged over samples, 



{W{C))t 



2m-l 



K„ 



(24) 



Furthermore, the free energy on the Nishimori line, apart 
from a nonsingular additive term, is equal to the Shan- 
non entropy of the distribution of magnetic monopoles, 
so that the latter is singular at p = p c . 



In principle, the RPGM could have what might be 
called a "gauge glass" phase. In this phase, the Wilson 
loop, averaged over thermal and quenched fluctuations, 
has area-law behavior, 



[(W(C))k] Kb ~exp[- K .Area(C)] 



(25) 



but the square of its thermal expectation value, averaged 
over quenched fluctuations, has perimeter-law behavior: 



[{W(C)) 2 K ] K ~ exp [-fj, • Perimeter(C*)] 



(26) 



This means that thermal fluctuations do not induce mag- 
netic disorder for each typical sample, but that the mag- 
netic fluctuations are large when we compare one sample 
to another. However, the identity eq. (24) shows that, 
along the Nishimori line K = K p , there can be no gauge 
glass phase. Since (W(C)) and (W(C)) 2 have the same 
average over samples, both order parameters cross from 
perimeter to area law at the same point on the Nishi- 
mori line. (Nishimori [10] used the analogous argument 
to show that there is no spin glass behavior in the RBIM 
along the Nishimori line.) 

Another useful identity that can be derived using the 
change of variable is 



\{W{C)) K ] = [(W(C)) K (W(C)) K: 



Since -1 < W(C) < 1, it follows that 



l(W(C)) K ] Kv < [ \{W(C)) Kp 



(27) 



(28) 



From this inequality, we may infer that if the point on the 
Nishimori line with concentration p is in the confinement 
phase, then the point (p, T) is in the confinement phase 
for any temperature T. (Again, the reasoning is exactly 
analogous to Nishimori's argument for the RBIM [10].) 
Since there is no gauge-glass behavior on the Nishimori 
line, if a point on the Nishimori line is in the confinement 
phase, then (W(C))k already exhibits area-law decay 
before averaging over samples. Therefore the right-hand 
side of eq. (28) shows area-law decay and so must the left- 
hand side. We conclude that, as for the RBIM, the phase 
boundary of the RPGM below the Nishimori line must 
either be vertical (parallel to the T axis as in Fig. 2a) or 
reentrant (tipping back toward the T axis as T decreases 
as in Fig. 2b). 



C. Further generalizations 

In d dimensions, the magnetic order parameter of 
the RBIM explores whether a thermally excited chain 
E' of codimension 1 (domain walls) deviates far from 
a quenched codimension- 1 chain E (antiferromagnctic 
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bonds), where both E and E' have the same codimension- 
2 boundary (the Ising vortices). Similarly, the RPGM 
can be defined in d dimensions, and its Wilson-loop order 
parameter probes whether a thermally excited chain E' of 
codimcnsion 2 (flux tubes) deviates far from a quenched 
codimension-2 chain E (Dirac strings), where both E and 
E' have the same codimension-3 boundary (the magnetic 
monopoles). 

This concept admits further generalizations. In d- 
dimensions, we may consider the lattice theory of a 
"rank-r antisymmetric tensor field" with quenched dis- 
order. Then variables reside on the r-cells of the lat- 
tice, and the Hamiltonian is expressed in terms of a field 
strength defined on (r + l)-cells. The sign of the coupling 
is determined by a random variable r taking values ±1 
on (r + l)-cclls; cells with the "wrong sign" have con- 
centration p. On the dual lattice, t corresponds to a 
codimcnsion- (r + 1) chain E, and the excited cells to 
a codimension-(r + 1) chain E', where E and E 1 are 
bounded by the same codimcnsion- (r + 2) chain of frus- 
tration. An operator analogous to the Wilson loop can 
be defined that detects the flux through the dimension- 
(r + 1) "surface" bounded by a dimension- r "loop" C; 
this operator serves as the order parameter for an order- 
disorder transition. The order parameter probes whether 
the thermally fluctuating codimension-(r+l) chain E' de- 
viates far from the quenched codimcnsion- (r + 1) chain 
E. 

For any d and r, the model has enhanced symmetry 
on the Nishimori line, where K = K p . Properties of the 
model on this line can be derived, analogous to those 
discussed above for the RBIM and the RPGM. 



III. ACCURACY THRESHOLD FOR QUANTUM 
MEMORY 

How the RBIM and RPGM relate to the performance 
of topological quantum memory was extensively dis- 
cussed in [9]. Here we will just briefly reprise the main 
ideas. 



A. Toric codes 

Quantum information can be protected from decoher- 
ence and other possible sources of error using quantum 
error-correcting codes [6,7] and fault-tolerant error re- 
covery protocols [8] . Topological codes (or surface codes) 
are designed so that the quantum processing needed to 
control errors has especially nice locality properties [4,5] . 

Specifically, consider a system of 2L 2 qubits (a qubit 
is a two-level quantum system), with each qubit residing 
at a link of an L x L square lattice drawn on a two- 
dimensional torus. (Other examples of surface codes, in- 
cluding codes defined on planar surfaces, are discussed 
in [9].) This system can encode two qubits of quantum 



information that are well protected from noise if the er- 
ror rate is low enough. The two-qubit code space, where 
the protected information resides, can be characterized 
as a simultaneous eigenspace with eigenvalue one of a 
set of check operators (or "stabilizer generators"); check 
operators are associated with each site and with each el- 
ementary cell (or "plaqucttc") of the lattice, as shown in 
Fig. 3. We use the notation 

'-(I !)■ x =(° °) ■ ,29 > 
?)•*-(; -°.) (»> 

for the 2x2 identity and Pauli matrices. The check oper- 
ator at site i acts nontrivially on the four links that meet 
at the site; it is the tensor product 

Xi = ®i 3s X t (31) 

acting on those four qubits, times the identity acting on 
the remaining qubits. The check operator at plaquette P 
acts nontrivially on the four links contained in the pla- 
quette, as the tensor product 

Zp = ®i e pZi , (32) 

times the identity on the remaining links. 





z 






z 




z 
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X 













FIG. 3. The check operators of the toric code. Each pla- 
quette operator is a tensor product of Z's acting on the four 
links contained in the plaquette. Each site operator is a ten- 
sor product of X's acting on the four links that meet at the 
site. 

The check operators can be simultaneously diagonal- 
ized, and the toric code is the space in which each check 
operator acts trivially. Because of the periodic bound- 
ary conditions on the torus, the product of all L 2 site 
operators or all L 2 plaquette operators is the identity 
- each link operator occurs twice in the product, and 
X 2 = Z 2 = I . There are no further relations among 
these operators; therefore, there are 2 • (L 2 — 1) inde- 
pendent check operators constraining the 2L 2 qubits in 
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the code block, and hence two encoded qubits (the code 
subspace is four dimensional). 

Since the check operators are spatially local, it is useful 
to think of a site or plaquette where the check operator 
has the eigenvalue —1 as the position of a localized exci- 
tation or "defect." The code space contains states with 
no defects, which are analogous to vacuum states of a Z 2 
gauge theory on the torus: Zp = 1 means that there is 
no Z2 magnetic flux at plaquette P, and = 1 means 
that there is no Z 2 electric charge at site i. (This Z 2 
gauge theory on the two-torus should not be confused 
with the three-dimensional Z2 gauge theory, described in 
Sec. IIIC, that arises in the analysis of the efficacy of 
error correction!) 
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FIG. 4. Site defects and plaquette defects in the toric code. 
Applied to the code space, Z's acting on a connected chain of 
links (darkly shaded) create site defects (electric charges) at 
the ends of the chain. Similarly, X's applied to a connected 
chain of dual links (lightly shaded) create plaquette defects 
(magnetic fluxes) at the ends of the chain. 

Consider applying to the vacuum state an operator 
that is a tensor product of Pauli matrices {Zg} acting 
on each of a set of links forming a connected chain {£}. 
This operator creates isolated site defects at the ends of 
the chain. Similarly, if we apply to the vacuum a tensor 
product of Pauli matrices {X^} acting on a connected 
chain of the dual lattice, isolated plaquette defects are 
created at the ends of the chain, as in Fig. 4. A general 
"Pauli operator" (tensor product of Pauli matrices) can 
be expressed as tensor product of Xg's and Ii's times a 
tensor products of Zgs and I/s; this operator preserves 
the code space if and only if the links acted upon by Z's 
comprise a cycle of the lattice (a chain with no bound- 
ary) and the links acted upon by X's comprise a cycle of 
the dual lattice. 

Cycles on the torus are of two types. A homologically 
trivial cycle is the boundary of a region that can be tiled 
by plaquettes. A product of Z's acting on the links of the 
cycle can be expressed as a product of the enclosed pla- 
quette operators, which acts trivially on the code space. 
A homologically nontrivial cycle wraps around the torus 



and is not the boundary of anything. A product of Z's 
acting on the links of the cycle preserves the code space, 
but acts nontrivially on the encoded quantum informa- 
tion. Associated with the two fundamental nontrivial cy- 
cles of the torus are encoded operations Z\ and Z 2 acting 
on the two encoded qubits. Similarly, associated with the 
two dual cycles of the dual lattice are the corresponding 
encoded operations X\ and X 2 , as shown in Fig. 5. 



(a) 



X 



2 
(b) 



FIG. 5. Basis for the operators that act on the two encoded 
qubits of the toric code, (a) The encoded Z\ is a tensor prod- 
uct of Z's acting on lattice links comprising a cycle of the 
torus, and the encoded X\ is a tensor product of X's acting 
on dual links comprising the complementary cycle, (b) Z 2 and 
X2 are defined similarly. 

A general error acting on the code block can be ex- 
panded in terms of Pauli operators. Therefore, we can 
characterize the efficacy of error correction by consider- 
ing how well we can protect the encoded state against 
Pauli operator errors. With the toric code, X errors (bit 
flips) and Z errors (phase flips) can be corrected inde- 
pendently; this suffices to protect against general Pauli 
errors, since a Y error is just a bit flip and a phase flip 
acting on the same qubit. We may therefore confine our 
attention to Z errors; the X errors may be dealt with in 
essentially the same way, but with the lattice replaced by 
its dual. 



B. Perfect measurements and the random-bond 
Ising model 

To be concrete, suppose that the Z errors are indepen- 
dently and identically distributed, occuring with proba- 
bility p on each qubit. Noise produces an error chain E, 
a set of qubits acted upon by Z. To diagnose the errors, 
the code's local check operators are measured at each 
lattice site, the measurement outcomes providing a "syn- 
drome" that we may use to diagnose errors. However, the 
syndrome is highly ambiguous. It does not completely 
characterize where the errors occured; rather it only in- 
dicates whether the number of damaged qubits adjacent 
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to each site is even or odd. That is, the measurement 
determines the boundary dE of the error chain E. 

To recover from the damage, we choose a recovery 
chain E' that has the same boundary as the measured 
boundary of E, and apply Z to all the qubits of E' . Since 
dE = dE', the chain D = E + E' is a cycle with dD = 0. 
Now, if D is homologically trivial, then the recovery suc- 
cessfully protects the two encoded qubits — the effect of 
the errors together with the recovery step is just to ap- 
ply a product of check operators, which has trivial action 
on the code space. But if D is homologically nontrivial, 
then recovery fails — the encoded quantum information 
suffers an error. 

Error recovery succeeds, then, if we can guess the ho- 
mology class of the randomly generated chain E, know- 
ing only its boundary dE — we succeed if our guess 
E' = E + D differs from E by a homologically trivial cy- 
cle D. If the error rate p is below a certain critical value 
p c called the accuracy threshold, it is possible to guess 
correctly, with a probability of failure that approaches 
zero for a sufficiently large linear size L of the lattice. 
But if p is above p c , the failure probability approaches 
a nonzero constant as L — > oo. The numerical value of 
p c is of considerable interest, since it characterizes how 
reliably quantum hardware must perform for a quantum 
memory to be robust. 

Let prob(-E) denote the probability that the error chain 
is E, and let prob[(E + D)\E] denote the normalized con- 
ditional probability for error chains E' = E+D that have 
the same boundary as E. Then, the probability of error 
per qubit lies below threshold if and only if, in the limit 
L — ► oo, 

^prob(-E)- prob[(£ + D)\E)= . (33) 

E D nontrivial 

Eq. (33) says that error chains that differ from the ac- 
tual error chain by a homologically nontrivial cycle have 
probability zero. Therefore, the outcome of the measure- 
ment of the check operators is sure to point to the correct 
homology class, in the limit of an arbitrarily large code 
block. 

This criterion is identical to the criterion for long-range 
order in the two-dimensional RBIM, along the Nishi- 
mori line. The error chain E can be identified with the 
chain of antifcrromagnetic bonds of a sample, bounded 
by Ising vortices that are pinned down by the measure- 
ment of the local check operators. The ensemble of all 
the chains {E'} with a specified boundary can be inter- 
preted as a thermal ensemble. If the temperature T and 
the error rate p obey Nishimori's relation, then the chain 
E' and the chain E have the same bond concentration. 
At low temperature along the Nishimori line, the cycle 
D = E + E' contains no large connected loops for typi- 
cal samples and typical thermal fluctuations — the spin 
system is magnetically ordered and error recovery suc- 
ceeds with high probability. But at higher temperature, 
the quenched chain E and the thermal chain E' fluctu- 
ate more vigorously. At the Nishimori point, D contains 



loops that "condense," disordering the spins and com- 
promising the effectiveness of error correction. Thus, 
the critical concentration p c at the Nishimori point of 
the two-dimensional RBIM coincides with the accuracy 
threshold for quantum memory using toric codes (where 
p c is the largest acceptable probability for either an X 
error or a Z error). 

The optimal recovery procedure is to choose a recov- 
ery chain E' that belongs to the most likely homology 
class, given the known boundary of the chain dE' = dE. 
For p < pc, the probability distribution has support on a 
single class in the limit L — > oo, and the optimal recov- 
ery procedure is sure to succeed. In the language of the 
RBIM, for a given sample with antiferromagnetic chain 
E, a chain E' of excited bonds can be classified according 
to the homology class to which the cycle D = E + E' be- 
longs, and a free energy can be defined for each homology 
class. For p < p c along the Nishimori line, the trivial ho- 
mology class has lowest free energy, and the free energy 
cost of choosing a different class diverges as L — ► oo. 

An alternative recovery procedure is to choose the sin- 
gle most likely recovery chain E' , rather than a chain 
that belongs to the most likely class. In the language of 
the RBIM, this most likely recovery chain E' for a given 
sample is the set of excited links that minimizes energy 
rather than free energy. This energy minimization proce- 
dure is sure to succeed if the error rate is p < p c o, where 
Pco is the critical bond concentration of the RBIM at 
T — 0. Since minimizing energy rather than free energy 
need not be optimal, we see that p c o < Pc- However, the 
energy minimization procedure has advantages: it can be 
carried out efficiently using the Edmonds perfect match- 
ing algorithm [19,20], and without any prior knowledge 
of the value of p. 



C. Imperfect measurement and the 
random-plaquette gauge model 

But the RBIM applies only to an unrealistic situation 
in which the outcomes of measurements of check oper- 
ators are known with perfect accuracy. Since these are 
four-qubit measurements, they must be carried out with 
a quantum computer and arc themselves degraded by 
noise. To obtain reliable information about the posi- 
tions of the Ising vortices, we must repeat the measure- 
ments many times, assembling a measurement history 
from which we can infer the world lines of the vortices in 
three-dimensional spacetime. 

To visualize the world lines in three dimensions, con- 
sider a three-dimensional simple cubic lattice on T 2 x R, 
where T 2 is the two-torus and R is the real line. The error 
operation acts at each integer-valued time t, and check 
operators are measured between each t and t+1. Qubits 
in the code block are associated with timelikc plaquettes, 
those lying in the tx and ty planes. A qubit error that oc- 
curs at time t is associated with a horizontal (spacelike) 
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bond that lies in the time slice labeled by t. An error in 
the measurement of a check operator at site j between 
time t and time t+ 1 is associated with the vertical (time- 
like) bond connecting site j at time t and site j at time 
t + 1. Qubit errors on horizontal bonds occur with prob- 
ability p, and measurement errors on vertical links occur 
with probability q. The set of all errors, both horizontal 
and vertical, defines a one-chain E, shown darkly shaded 
in Fig. 6. The set of all syndrome measurements with 
nontrivial outcomes (those where the observed value of 
the check operator is —1 rather than +1) defines a (ver- 
tical) one-chain S, shown lightly shaded in Fig. 6. The 
chains E and S share the same boundary; therefore the 
(possibly faulty) measurements of the check operators re- 
veal the boundary of the error chain E. 




FIG. 6. An error history shown together with the syndrome 
history that it generates, for the toric code. For clarity, the 
three-dimensional history of the two-dimensional code block 
has been compressed to two dimensions. Qubits reside on 
plaquettes, and four-qubit check operators are measured at 
each vertical link. Links where errors have occured are darkly 
shaded, and links where the syndrome is nontrivial are lightly 
shaded. Errors on horizontal links indicate where a qubit 
flipped between successive syndrome measurements, and er- 
rors on vertical links indicate where the syndrome measure- 
ment was wrong. Vertical links that are shaded both lightly 
and darkly are locations where a nontrivial syndrome was 
found erroneously. The chain S of lightly shaded links (the 
syndrome) and the chain E of darkly shaded links (the errors) 
both have the same boundary. 

Error recovery succeeds if we can guess the homology 
class of the error chain E, given knowledge of its bound- 
ary dE; that is, we succeed if our guess E' = E + D 
differs from E by a cycle D that is homologically trivial 
on T 2 x R. Thus, the accuracy threshold can be mapped 
to the confinement-Higgs transition of the RPGM. The 
error one-chain E on the dual lattice becomes the set of 
wrong-sign plaquettes on the lattice; its boundary points 
are magnetic monopoles, whose locations are determined 
by the measurements of local check operators. Since q 
need not equal p, the gauge model can be anisotropic 



- on the original lattice, the concentration of space- 
like wrong-sign plaquettes is q (spacelike plaquettes are 
dual to timelike bonds) and the concentration of timclikc 
wrong-sign plaquettes is p (timelike plaquettes are dual to 
spacelikc bonds). The ensemble of error chains {E'} that 
have the same boundary as E becomes the thermal en- 
semble determined by an anisotropic Hamiltonian, with 
the coupling -ftT spaco on spacelike plaquettes obeying the 
Nishimori relation ifspaco = K q and the coupling -Ktime 
on timelikc plaquettes the relation if time = K p . 

For small p and q, the cycle D = E + E' contains 
no large connected loops for typical samples and typical 
thermal fluctuations — the gauge system is magnetically 
ordered and error recovery succeeds with high probabil- 
ity. But there is a critical curve in the (p, q) plane where 
the magnetic flux tubes "condense," magnetically disor- 
dering the system and compromising the effectiveness of 
error correction. For the sort of error model described in 
[9], the qubit error rate and the measurement error rate 
are comparable, so the isotropic model with p = q pro- 
vides useful guidance. For that case, the critical concen- 
tration p c at the Nishimori point of the three-dimensional 
RPGM coincides with the accuracy threshold for quan- 
tum memory using toric codes (where p c is the largest 
acceptable probability for an X error, a Z error, or a 
measurement error). In the extreme anisotropic limit 
q — ► 0, flux on spacelike plaquettes is highly suppressed, 
and the timelike plaquettes on each time slice decouple, 
with each slice described by the RBIM. 

For both the 2D RBIM and the 3D (isotropic) RPGM, 
we may infer (as Nishimori argued for the RBIM [10]) 
that the phase boundary lies in the region p < p c , i.e., 
does not extend to the right of the Nishimori point. From 
the perspective of the error recovery procedure, this prop- 
erty reflects that the best hypothesis about the error 
chain, when its boundary is known, is obtained by sam- 
pling the distribution prob[(_E + D)\E\. Thus, for each 
value of p, the fluctuations of D are best controlled (the 
spins or gauge variables are least disordered) by choosing 
the temperature on the Nishimori line. For p > p c the 
magnetization of the 2D RBIM vanishes on the Nishi- 
mori line, and so must vanish for all T . A similar remark 
applies to the Wilson-loop order parameter of the 3D 
RPGM. 



In particular, the critical value of p on the T = axis 
(denoted p c0 ) provides a lower bound on p c . Rigorous 
arguments in [9] established that p c o > .0373 in the 2D 
RBIM and p c0 > .0114 in the 3D RPGM. (A similar 
lower bound for the 2D RBIM was derived by Horiguchi 
and Morita many years ago [21].) We have estimated the 
value of Pco using numerical simulations that we will now 
describe. 
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IV. NUMERICS 



A. Method 



For the RBIM in two dimensions (but not in higher di- 
mensions), and for the RPGM in three dimensions (but 
not in higher dimensions), it is numerically tractable to 
study the phase transition on the T = axis. Specifically, 
for the RBIM, we proceed as follows: Consider anLxi 
lattice on the torus, and generate a sample by choosing 



a random at each bond (where Tj. 



-1 occurs with 



probability p) . Consider, for this sample, the one-chain E 
on the dual lattice containing bonds with — — 1, and 
compute its boundary dE to locate the Ising vortices. 

Then, to find the ground state of the Hamiltonian for 
this sample, construct the one-chain E' of the dual lat- 
tice, bounded by the Ising vortices, with the minimal 
number of bonds. This minimization can be carried out 
in a time polynomial in L using the Edmonds perfect 
matching algorithm [19,20]. (If the ground state is not 
unique, choose a ground state at random.) Now exam- 
ine the one-cycle D = E + E' on the torus and compute 
whether its homology class is trivial. If so, we declare 
the sample a "success;" otherwise the sample is a "fail- 
ure." Repeat for many randomly generated samples, to 
estimate the probability of failure Pf a ii(p). 

We expect Pf a ii(p) to be discontinuous at p = p c o in 
the infinite volume limit. For p < p c g, large loops in D 
are heavily suppressed, so that Pf a i] falls exponentially to 
zero for L sufficiently large compared to the correlation 
length £. But for p > p c0 , arbitrarily large loops are not 
suppressed, so we anticipate that the homology class is 
random. Since there are four possible classes, we expect 
Pfaii to approach 3/4 as L — > oo. 

This expectation suggests a finite-size scaling ansatz 
for the failure probability. Let the critical exponent v$ 
characterize the divergence of the correlation length £ at 
the critical point p = p c o'. 



£ ~ |p-Pco| 



(34) 



For a sufficiently large linear size L of the sample, the 
failure probability should be controlled by the ratio L/£; 
that is, it is a function of the scaling variable 

x={p- Pco)L 1/,Jo . (35) 

Thus the appropriate ansatz is 

Pfaii ~ \f{x) , (36) 
where the function / has the properties 

lim f{x) = , lim f(x) = 1 . (37) 

x^ — oo x^oo 

Though the scaling ansatz should apply asymptotically 
in the limit of large L, there are systematic corrections 
for finite L that are not easily estimated. 



According to eq. (36), the failure probability at p = p c o 
has a universal value (3/4)/(0) that does not depend on 
L. Thus, by plotting Pf a ;i vs. p for various values of L, 
we can estimate p c o by identifying the value of p where 
all the curves cross. To find z/q, we observe that 



log 



fdP,, 



ail 



1 



log L + constant . 



(38) 



Hence, if we estimate the slope of Pf a ;i at p = p c o, we 
can extract vq from a linear fit to a plot of log(slope) vs. 
logL. 

The three-dimensional RPGM can be analyzed by the 
same method. A sample is generated by randomly choos- 
ing Tp on each plaquette of an L 3 cubic lattice on the 3- 
torus. The wrong-sign plaquettes define a one-chain E on 
the dual lattice, whose boundary defines the locations of 
the magnetic monopoles. The ground state of the sample 
is constructed by finding the one-chain E' with the same 
boundary that has the minimal length, and the one-cycle 
D = E + E 1 is examined to determine if it is homolog- 
ically trivial. Since there are eight homology classes on 
the 3-torus, the scaling ansatz becomes 
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Pfaii ~ o/(ar) , 



(39) 



and p c0 and vq are estimated as described above. 

For the RBIM in three dimensions, or the RPGM in 
four dimensions, E and E' become two-chains. To con- 
struct the ground state, then, we must find the minimal 
two-dimensional surface that has a specified boundary. 
Unfortunately, this problem is known to be NP-hard [22] 
and so appears to be computationally intractable. 

Detailed numerical studies of the two-dimensional 
RBIM in the vicinity of the Nishimori point have been 
done earlier by other authors [11,12], using methods that 
are not very effective at low temperature. The T = 
phase transition has been studied using methods related 
to ours [20,23], but with less numerical accuracy. As far 
as we know, numerical studies of the RPGM have not 
been previously attempted. 



B. Random-bond Ising model 

We measured Pf a ii by generating 10 6 samples for each 
value of L from 2 to 36, and for each value of p increas- 
ing in increments of .001 from .100 to .107; in addi- 
tion we generated 10 6 samples at L = 37,38,40,42 for 
p = .102, .103, .104. Values of P fai i for even L lie slightly 
but systematically above the values for odd L at the same 
p; therefore we analyzed the data for even and odd L sep- 
arately. Data for L = 16,20,24,28,32,36 are shown in 
Fig. 7, and data for L = 15,19,23,27,31,35 are shown 
in Fig. 8. Crudely, the point of concordance of the data 
sets provides an estimate of p c o, while the trend of the 
data with L determines the exponent fo- 
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FIG. 7. The failure probability Pf ai i as a function of the er- 
ror probability p for linear size L — 16, 20, 24, 28, 32, 36, in the 
two-dimensional random-bond Ising model. Each data point 
was generated by averaging 10 6 samples. 
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FIG. 8. The failure probability Pf ai i as a function of the er- 
ror probability p for linear size L — 15, 19, 23, 27, 31, 35, in the 
two-dimensional random-bond Ising model. Each data point 
was generated by averaging 10 6 samples. 

We did a global fit of the data to the form 

P fai] = A + Bx + Cx 2 , (40) 

where x = (p — p^L 1 /" , adopting a quadratic approx- 
imation to the scaling function f(x) in the vicinity of 
x = 0. (In the range of x we considered, the quadratic 
term is small but not quite negligible.) For even L rang- 
ing from 22 to 42, our fit found 



p c0 = .10330 ± .00002 
V{) = 1.49 ± .02 , 



(41) 



where the quoted errors are one-sigma statistical errors. 
For odd L ranging from 21 to 37, our fit found 



p c0 = .10261 ±.00003 
v = 1.46 ± .02 . 



(42) 



The discrepancy between the values of p c o for even and 
odd L indicates a nonncgligible finite-size effect. 



-0.01 0.01 0.02 

rescaled error rate 



FIG. 9. The failure probability Pfaii, with the nonuniversal 
correction of Eq. (43) subtracted away, as a function of the 
scaling variable x = (p — p c o)L 1 / l/ ° for the two-dimensional 
random-bond Ising model, where p c o and vq are determined 
by the best fit to the data. A two-sigma error bar is shown 
for each point. The data for values of L from 2 to 42 lie on 
a single line, indicating that the (small) scaling violations are 
well accounted for by our ansatz. 

On closer examination, we see evidence for small but 
detectable violations of our scaling ansatz in both the 
even and odd data sets. These violations are very well 
accounted for by the modified ansatz 



ail 



Cx 2 

even ^ 



Bx 
D, 

'odd 



(L even) , 
(L odd) , 



(43) 



which includes a nonuniversal additive correction to Pf a n 
at criticality, different for even and odd sizes. Fitting the 
modified ansatz to the data for even L ranging from 2 to 
42, we find 



Pco = .10309 ± .00003 , 
v = 1.461 ±.008, 
D cvon = 0.165 ± .002 , fi cvcn = 0.71 ± .01 . 



(44) 



Fitting to the data for odd L ranging from 3 to 37, we 
find 



p cQ = .10306 ± .00008 , 
vq = 1.463 ±.006, 
AxM = -.053 ±.003 , Modd = 2.1 ± .3 . 



(45) 



In Fig. (9) we show the data for all values of L and p; 
using the values of p c0} D, and \i found in our fits, 
we have plotted Pfaii, with the nonuniversal correction 
of Eq. (43) subtracted away, as a function of the scaling 
variable x = (p — p^L 1 ^ . All of the data lie on a single 
line, indicating that residual scaling violations are quite 
small. Furthermore, the agreement between the values 
of Pco and vq extracted from the even and odd data sets, 
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which were fit independently, indicates that our extrap- 
olation to large L is reasonable, and that the statistical 
errors in Eq. (44,45) do not seriously underestimate the 
actual errors in our measurement. A plausible conclusion 
is that 



p c0 = .1031 ± .0001 
vq = 1.46 ± .01 . 



(46) 



An earlier measurement reported by Kawashima and 
Rieger found [23] 



p c0 = .104±.001 , 
v = 1.30 ± .02 ; 



(47) 



their value of p c o, but not of vq, is compatible with ours. 
An important reason why our value of p c o has a smaller 
statistical error than theirs is that they computed a dif- 
ferent observable (the domain wall energy) for which the 
finite-size scaling analysis is more delicate than for the 
failure probability (another critical scaling exponent is 
involved) . 

In a recent study of the Nishimori point, Merz and 
Chalker found [12] 



p c = .1093 ± .0002 , 
v = 1.50 ± .03 . 



(48) 



L = 10 
L = 12 
L = 14 
L = 16 



0.029 0.0 

error rate p 



FIG. 10. The failure probability Pf a u as a function of the 
error probability p for linear size L = 10, 12, 14, 16, in the 
three-dimensional random-plaquette gauge model. Each data 
point was generated by averaging 10 6 samples. 



We did a global fit of the data to the form 



ffail = A + Bx + Cx 2 



(50) 



where x = (p — p c ^)L x l V( \ adopting a quadratic approx- 
imation to the scaling function f(x) in the vicinity of 
x = 0. For L ranging from 9 to 16, our fit found 



There is a clear discrepancy between the values of p c and 
Pco, in disagreement with the conjecture of Nishimori [13] 
and Kitatani [14]. Evidence for a reentrant phase dia- 
gram has also been found by Nobre [24], who reported 

p c0 = .1049 ± .0003 . (49) 

In principle, the phase transitions at T = and at the 
Nishimori point could be in different universality classes, 
so that the critical exponents vq and v could have differ- 
ent values. However, our measurement of v at T = 
is consistent with the value of v at the Nishimori point 
reported by Merz and Chalker [12]. 



C. Random-plaquette gauge model 

We measured Pf a n by generating 10 6 samples for each 
value of L from 9 to 14, and for each value of p in- 
creasing in increments of .0004 from .02805 to .03005; 
in addition we generated 10 6 samples at L = 15, 16 for 
p = .02845, .02925, .03005. Values of P iai i for even L lie 
slightly but systematically above the values for odd L at 
the same p; therefore we analyzed the data for even and 
odd L separately. Data for even L are shown in Fig. 10. 
Crudely, the point of concordance of the data sets pro- 
vides an estimate of p c o, while the trend of the data with 
L determines the exponent vq. 



p c0 = .02937 ± .00002 , u Q = 0.974 ± .026 (L even) , 
p c0 = .02900 ± .00001 , v Q = 1.025 ± .016 (L odd) , (51) 

where the quoted errors are one-sigma statistical errors. 
The results for even and odd L are incompatible, indi- 
cating a nonnegligible finite-size effect. 

0-22 I — i 1 1 1 1 1 i — I 




0.! L_l 1 1 1 1 1 L_ 

-0.02 -0.015 -0.01 -0.005 0.005 0.01 

rescaled error rate 



FIG. 11. The failure probability Pf a n as a function of the 
scaling variable x = (p — p^L 1 ^" for the random-plaquette 
gauge model, where p c o and vo are determined by the best fit 
to the data. A two-sigma error bar is shown for each point. 
The data for all even values of L from 10 to 16 lie on a single 
curve, indicating that scaling violations are small. 
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We believe that our analysis for even L is likely to be 
more reliable; finite size effects are enhanced for odd L, 
the case in which the failure probability is smaller. All of 
the even-L data are shown in Fig. 11, with Pf a n plotted 
as a function of x = (p — p c o)L 1 / u ° , where p c o and uq 
are determined by our fit. The data fit a single curve, 
indicating that scaling violations are small. (Scaling vi- 
olations are more discernable in the odd-L data set.) A 
reasonable conclusion is that 



PcO 



.0293 ± .0002 
1.00 ± .05 . 



(52) 



D. The failure probability at finite temperature 

Our numerical studies of the RBIM and the RPGM 
were restricted to the T = axis. We calculated the fail- 
ure probability to estimate the critical disorder strength 
Pco and the critical exponent vq. Here we will describe 
how the calculation of the failure probability could be 
extended to nonzero temperature. 

To calculate the failure probability in the zero- 
temperature RBIM, we generate a sample by specifying a 
one-chain E of antiferromagnetic links, and then we con- 
struct the one-chain E' of minimal energy with the same 
boundary as E. Failure occurs if the cycle D = E + E' 
is homologically nontrivial. 

At nonzero temperature we should construct E' to be- 
long to the homology class that minimizes free energy 
rather than energy. For a given sample with antiferro- 
magnetic one-chain E, the free energy F(E, h) of ho- 
mology class h is found by summing over domain wall 
one-chains {£"} such that E + E' e h: 



exp{-PF(E,h)] = Z(E,h)= Yl ' 

E'-.E+E'eh 



-0H E 



(53) 



where He denotes the Hamiltonian eq. (1) with antifer- 
romagnetic chain E. If the trivial homology class h = e 
has the lowest free energy, then the sample is a "success;" 
otherwise it is a "failure." We can estimate the failure 
probability -Pf a ii(p, T) by randomly generating many sam- 
ples, and determining for each whether it is a success or 
a failure. 

For the random bond Ising model on a torus, the sum 
eq. (3) includes only the chains E' such that E + E' is in 
the trivial homology class. To sum over the class h, we 
can augment E by adding to it a representative of h. For 
each h, we can compute 



Z(E,h) 
~ZjE^e) 



= exp[-f3(F(E,h)-F(E,e))] 



(54) 



the sample E is a success if this ratio of partition func- 
tions is less than one for each h ^ e. 

The ratio is the thermal expectation value (Oh)k of 
an observable Oh that "inserts a domain wall" wrapping 



around a cycle C representing h. That is, the effect of 
Oh is to flip the sign of the bond variable for each 
bond (ij) in C: 



O h = exp 



-IK ]T • 



ij Si Sj 



(55) 



In principle, we could measure (Oh)K by the Monte Carlo 
method, generating typical configurations in the thermal 
ensemble of He, and evaluating Oh in these configura- 
tions. Unfortunately, this method might not produce an 
accurate measurement, because the configurations that 
dominate (Oh)K may be exponentially rare in the ther- 
mal ensemble — a configuration with excited bonds on 
C can have an exponentially large value of Oh that over- 
comes exponential Boltzmann suppression. 

One solution to this problem is to express 
Z(E,h)/Z(E,e) as a product of quantities, each of 
which can be evaluated accurately by Monte Carlo. Let 
{e = P ,Pi,P 2 , . . . Pk-i,Pk = C} be a sequence of open 
chains interpolating between the empty chain and the 
cycle C, where Pj+i — Pj contains just a single bond. 
We may write 



Z{E,h) Z(E,P X ) Z(E,P 2 



Z{E,e) Z(E,P ) Z{E,P X ) 



Z{E,P k ) 
Z(E,P k ^) 



(56) 



Each ratio Z(E, Pj + i)/Z(E, Pj) is the expectation value 
of an operator that acts on a single bond, evaluated in 
the thermal ensemble of the Hamiltonian with antiferro- 
magnetic bonds on the chain E + Pj; this expectation 
value can be evaluated by Monte Carlo with reasonable 
computational resources. (For an application of this trick 
in a related setting, see [25].) 

Using this method, we can determine whether 
Z(E,h)/Z(E,e) exceeds one for any h ^ e and hence 
whether the sample E is a success or a failure. Gen- 
erating many samples, we can estimate Pf a n(p,T). In 
principle, then we can calculate the failure probability 
for the optimal recovery scheme, in which p and T obey 
Nishimori's relation. By a similar method, we can calcu- 
late the failure probability for the RPGM. However, we 
have not attempted this calculation. 



V. CONCLUSIONS 

The three-dimensional random-plaquette gauge model, 
and the analogous antisymmetric tensor models in higher 
dimensions, provide new examples of multicritical points 
with strong disorder. These models have phase diagrams 
that qualitatively resemble the phase diagram of the two- 
dimensional random-bond Ising model. 

Our results indicate that the boundary between the 
ferromagnetic and paramagnetic phases of the RBIM is 
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reentrant rather than vertical below the Nishimori line. 
If the disorder strength p satisfies p c o < p < p c , then the 
ground state of the spin system does not have long-range 
order. As the temperature T increases with p fixed, long- 
range order is first restored, and then lost again as the 
temperature increases further. At T = the spins arc 
frozen in a disordered state driven by quenched random- 
ness. But apparently this ground state is entropically 
unfavorable — at low but nonzero temperature typical 
states in the thermal ensemble have long-range ferromag- 
netic order. 

This behavior seems less remarkable when considered 
from the viewpoint of our error recovery protocol. For 
given p and a specified error syndrome, the recovery 
method with optimal success probability proceeds by in- 
ferring the most likely homology class of errors consistent 
with the syndrome. There is no a priori reason for the 
most likely single error pattern (the ground state) to be- 
long to the most likely error homology class (the class 
with minimal free energy) even in the limit of a large 
sample. Our numerical results indicate that for error 
probability p such that p c0 < p < p c , the probability that 
the ground state does not lie in the most likely homology 
class remains bounded away from zero as L — > oo. 

In our numerical studies of the RBIM and RPGM at 
zero temperature, we have computed a homological ob- 
servable, the failure probability. This observable has 
advantages over, say, the domain wall energy, because 
it obeys a particularly simple finite-size-scaling ansatz. 
Therefore, we have been able to determine the critical 
disorder strength p C Q and the critical exponent i>$ to good 
accuracy with relatively modest computational resources. 

Not surprisingly, our numerical values for p c0 are no- 
tably larger than rigorous lower bounds derived using 
crude combinatoric arguments in [9]: p c o « .1031 com- 
pared with the bound p c a > .0373 in the RBIM, and 
p c0 w .0293 compared with p c0 > .0114 in the RPGM. 

The zero-temperature critical disorder strength p c0 is a 
lower bound on the value of the critical disorder strength 
p c along the Nishimori line, and of special interest be- 
cause of its connection with the accuracy threshold for ro- 
bust storage of quantum information. Our result means 
that stored quantum data can be preserved with arbi- 
trarily good fidelity if, in each round of syndrome mea- 
surement, qubit errors and syndrome measurement er- 
rors are independently and identically distributed, with 
error probability per qubit and per syndrome bit both 
below 2.9%. For qubit errors and measurement errors 
occuring at differing rates, an accuracy threshold could 
be inferred by analyzing an anisotropic random-plaquette 
gauge model, with differing disorder strength for horizon- 
tal and vertical plaquettes. Relating these threshold error 
rates to fidelity requirements for quantum gates requires 
further analysis of the sort discussed in [9] . 

We have also measured the critical exponent v§ that 
controls the divergence of the correlation length as p ap- 
proaches pco, finding v ss 1.46 in the RBIM and v s=s 1.0 
in the RPGM. The value of v is also relevant to the 



efficacy of quantum error correction — through its con- 
nection with finite-size scaling, v$ determines how large 
the code block should be to achieve a specified storage 
fidelity, for p less than but close to p C Q. 

Quantum computers are believed to be more power- 
ful than classical computers — classical computers are 
unable to simulate quantum computers efficiently The 
accuracy threshold for quantum memory is a fascinating 
phase transition, separating a low-noise quantum phase 
from a high-noise classical phase. In this paper, we have 
described one way to analyze this phase transition us- 
ing methods from traditional statistical physics. Further- 
more, the connection with quantum memory provides an 
enlightening new perspective on local spin and gauge sys- 
tems with strong quenched disorder. 
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